Fast decay of the velocity autocorrelation function in dense shear 
flow of inelastic hard spheres 
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Abstract. - We find in complementary experiments and event driven simulations of sheared 
inelastic hard spheres that the velocity autocorrelation function decays much faster than 
t'^^"^ obtained for a fluid of elastic spheres at equilibrium. Particle displacements are measured 
in experiments inside a gravity driven flow sheared by a rough wall. The average packing fraction 
obtained in the experiments is 0.59, and the packing fraction in the simulations is varied between 
0.5 and 0.59. The motion is observed to be diffusive over long times except in experiments where 
there is layering of particles parallel to boundaries, and diffusion is inhibited between layers. 
Regardless, a rapid decay of ijjit) is observed, indicating that this is a feature of the sheared 
dissipative fluid, and is independent of the details of the relative particle arrangements. An 
important implication of our study is that the non-analytic contribution to the shear stress may 
not be present in a sheared inelastic fluid, leading to a wider range of applicability of kinetic 
theory approaches to dense granular matter. 



A significant feature of the dynamics of dense fluids is 
the effect of correlations on transport coefficients [1,2]. 
For a gas of elastic particles in the dilute limit, the trans- 
port coefficients are calculated using the molecular chaos 
approximation [3] , that the two-particle velocity distribu- 
tion is the product of the single particle velocity distribu- 
tion functions. It is known that this procedure cannot be 
extended to higher densities due to the effect of correla- 
tions. For a D dimensional elastic fluid, the velocity auto- 
correlation function iplt) decays as t~^l'^ , and this power 
law decay is referred to as the "long time tail" in the ve- 
locity autocorrelation function [4]. This was recognized 
as a consequence of the diffusive transport of momentum 
proportional to the square of the wave vector in a sys- 
tem in which momentum and energy are conserved. This 
leads to a divergence of the viscosity in two dimensions. 
In three dimensions, there is a non-analytic contribution 
to the shear stress proportional to 7'^/^ in limit of zero 
shear rate [1]; this correction is larger than the Burnett 
terms (proportional to 7^ in the constitutive relation for 
the shear stress obtained by the Chapman-Enskog proce- 
dure. These correlations have been anticipated to impact 
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extention of kinetic theory to dense granular flows since 
its first application to dilute granular gases, and have been 
also seen in uniform granular flows [5]. Because dense 
flows are more prevalent, such observations have led to 
the belief that kinetic theory approaches are severely lim- 
ited. 

However, more recently, it has been speculated that 
the decay of i\){t) in sheared granular flows could be 
faster [6,7]. For a fluid under shear, a steady state can be 
obtained only when there is an energy dissipation mech- 
anism, since there is a continuous increase in the fluctu- 
ating (thermal) energy of the fluid particles, and the rate 
of increase is equal to the product of the stress and the 
strain rate. The two most widely studied mechanisms are 
thermostats (where there is a drag force on the particles) 
and inelastic collisions. In a sheared elastic fluid in the 
absence of thermostats, there is a continuous increase in 
the thermal energy (temperature) of the particles. Due to 
this, it is not possible to calculate the decay of velocity 
autocorrelation functions at steady state; the decay of the 
autocorrelation function will depend both on the initial 
temperature and the rate of heating, and it will not be 
invariant with respect to translation in time. The steady 
state is completely defined only if we specify both the en- 
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ergy production mechanism (mean shear) and the energy 
dissipation mechanism (inelastic coUisions or thermostat) . 

The velocity autocorrelation function in a thermostated 
system with inelastic collisions and random noise has been 
examined by Puglisi et al [9], and they have found that 
the decay is not exponential, but rather has a behaviour 
similar to a stretched exponential. This is in contrast to 
the algebraic decay in an elastic fluid at equilibrium in the 
absence of shear [4]. The reason for this difference could 
be the following. In an elastic fluid, the long time tails 
in the velocity correlation functions is due to the diffusive 
transport of the transverse momentum in the fluid [1]. 
In a thermostated system, momentum is not a conserved 
variable locally, since there are drag and fluctuating forces 
acting on the individual particles. This could result in a 
faster decay in the velocity autocorrelation functions. 

Here, we analyse a class of sheared dissipative fluids 
with inelastic collisions, where momentum is conserved in 
collisions, and energy is dissipated. Due to this, energy 
is now a non-conserved variable, and we examine the ef- 
fect of mean shear and energy dissipation on the velocity 
autocorrelation function. It has been suggested [6,7] that 
the velocity autocorrelation function in a sheared inelastic 
fluid should decay much faster than the t~^l'^ power law 
decay in an elastic fluid at equilibrium. 

Kumaran [7] incorporated correlations using the ring- 
kinetic equation for a dense gas to obtain a faster decay 
proportional to in three dimensions. A later more 

detailed calculation [8], which accounted for the rotation 
of the wave vector with the mean shear for calculating the 
wave-number cut-offs for perturbations modulated in the 
flow direction, indicated that the decay rate of the velocity 
autocorrelation function should be anisotropic, and should 
scale as t~^^l^ in the flow plane, and as in the vortic- 
ity direction. These results are in contrast with the power 
law decay obtained by Otsuki and Hayakawa [6] in 

three dimensions using mode coupling theory. Either of 
these corrections would render the viscosity finite in two 
dimensions, and there is no non-analytic contribution to 
the viscosity proportional to 7'^/^ in three dimensions, sug- 
gesting that calculations based on kinetic theory would 
have a greater range of validity, at least upto Burnett or- 
der in both two and three dimensions, for granular flows 
in comparison to elastic fluids. 

We investigate the velocity autocorrelation function in a 
sheared granular flow in complementary experiments and 
simulations. In experiments, we produce a sheared gran- 
ular flow with a rough boundary, and measure particle 
motion in three dimensions away from side walls using an 
interstitial liquid index matching technique. Further, to 
decouple the effect of the boundary on the correlations, to 
examine effect of volume fraction, and to examine a sim- 
pler granular system undergoing linear shear, we perform 
event driven (ED) simulations of inelastic spheres with 
Lees-Edwards boundary conditions. Therefore, in the ex- 
perimental system the flow occurs past boundaries which 
are held fixed and the granular volume fraction can vary 
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Fig. 1: (a) Schematic diagram of the granular flow cell used 
in the experiments. A layer of glass beads is glued to the 
right side wall in order to roughen it. The observation area 
is indicated by the dashed box. (b) The mean velocity and 
the strain rate 7 across the shear plane, (c) The number of 
particles per unit volume across the shear plane shows layering 
near the side walls. 



within the flow. In contrast, the simulations are performed 
with fixed volume fraction conditions. Independent of 
these details, the experiments and simulations show rapid 
decay of autocorrelations in a dense packing regime with 
volume fraction as high as 0.59. Further, the particles 
in both experiments and simulations do not show order- 
ing in the plane perpendicular to the direction of shear, 
as evidenced by the in-plane hexagonal order parameter. 
Layering occurs near the boundaries in the experiments, 
while there is no order even in the gradient direction in 
the simulations because we use periodic boundary condi- 
tions. Despite these differences, ■0(t) shows a decay much 
faster than t'^l"^ in both cases. At the onset, we note that 
because of considerable challenges in measuring and track- 
ing experimental particles in the bulk, the exact exponent 
is beyond our capability as it requires tracking particles 
longer than they are in our system. 

The experimental system consists of a vertical glass pipe 
with a rectangular cross section filled with glass beads with 
diameter d — 1.0 ± 0.1 mm, and density pg = 2.5 x 10'^ kg 
m^"^ driven by gravity. A layer of glass beads is glued to 
one of the sides to roughen it as shown in Fig. 1(a). The 
interstitial space between the grains is filled with a liquid 
with the same refractive index as the glass beads [10, 11]. 
The liquid has density pi = 1.0 x 10'^ kg m~'^, and viscosity 
v = 2.2 X lO^'^ kg m~^s. A dye added to the liquid is illu- 
minated by a light sheet and imaged from an orthogonal 
direction. Thus the particles appear dark against a bright 
background, and a centroid algorithm is used to find the 
position of the particles to within a twentieth of a particle 
diameter. By switching the light and image planes, we 
can measure the particle motion in the shear and in the 
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Fig. 2: (a) Velocity auto-correlations function versus time 
scaled with strain rate obtained from the experiments, {i — 
X — O, y — o, z — <>). tpii is observed to become more isotropic 
(right to left) as one moves away from the shearing wall where 
particles form layers. The four layers are centered around 
1.57d, 2.54d, 3.49d and 4.48d from the rough wall. V(i) mea- 
sured in the experiment plotted in linear scale (b) and logscale 
(c). ^{t) decay faster in sheared flow (open symbol) compared 
with unsheared flows at various flow rates in Ref. [5] (flUed 
symbol). Lines with slope —5/2 (dash dot), —7/2 (dot) and 
— 15/4 (dash) are also added to guide the eye. An exponential 
fit (solid line) is also added for comparision. 



vorticity planes. 

Figure 1(b) shows the measured mean flow velocity 
< Vx > &s a function of horizontal distance y in the shear 
plane. To minimize the effect of the front and back walls, 
the data is obtained for 20d < z < 40d. The rough bound- 
ary is observed to slow the flow compared to the smooth 
boundary. At these measured velocities, the Reynolds 
number Re = Vxd/v is about 10~^, and the ratio of the 
viscous drag of the grain to the gravitational force is es- 
timated to be less than 10~^. The strain rate 7 given by 
the gradient of the mean velocity is shown in Fig. 1(b) 
as well. The maximum strain rate 0.25s~^, is compara- 
ble in magnitude to a previous granular shear study [11] 
where the drag friction experienced by the boundary was 
measurably unaffected by the presence of the same inter- 
stitial liquid. Thus we conclude that the lubrication forces 



between grains do not affect the grain motion in the re- 
gion where we discuss the correlation properties. The side 
walls also impose layering as can be observed from the os- 
cillations in the plot of number of particles in a volume 
plotted across the horizontal crossection in Fig. 1(c). 
The oscillations decrease with distance from the bound- 
aries, and can be used to measure any influence of grain 
order on correlations. By taking the average of n^d~^ and 
multiplying by the volume of sphere, we obtain that the 
average volume fraction of the grains 4> in the observation 
window is 0.59. 



ED simulations are performed in a three dimensional 
box with 500 smooth inelastic spheres of diameter d sub- 
jected to a rate of deformation field G = '^e^^y. The co- 
ordinate system is similar to that in the experiments, and 
the flow is in the x direction, the velocity gradient in the y 
direction and the z co-ordinate is perpendicular to the flow 
plane. The box size is adjusted to obtain the desired vol- 
ume fraction. Simulations in three dimensions at large vol- 
ume fractions (upto 0.6) and with coefficient of restitution 
as low as 0.7 have been carried out using the ED method 
by several authors [12]. In the ED procedure, particle tra- 
jectories are extrapolated forward in time and the earliest 
possible collision is identified. All particle positions are 
advanced by this time interval, and the post-coUisional 
velocity of the colliding particles is updated along the line 
joining centers of the particles with a factor — e„ times 
the pre-collisional relative velocity, where e„ is the coef- 
ficient of restitution. The simulations were initiated with 
the particles arranged in an FCC lattice, and with veloc- 
ities chosen from a random distribution with variance 1. 
The simulation was first allowed to proceed for 2 x 10^ 
collisions per particle in order to reach a steady state, 
and the averaging was carried out over another 2 x 10^ 
collisions per particle. Shear flow is imposed using Lees- 
Edwards periodic boundary conditions. More details on 
the simulation methods can be found in Ref. [13]. For 
an inelastic fluid under shear, the strain rate 7 and the 
'granular temperature' T (mean square of the fluctuating 
velocities of the particles) are related through the energy 
balance condition, that the rate of shear production and 
inelastic dissipation of energy of the particles are equal 
at steady state. We carry out simulations with smooth 
spheres with e„ = 0.9, and with the volume fraction up 
to (0 = 0.59). For these parameters, we have taken care 
to ensure that there is no inelastic collapse [14], and there 
are no particle overlaps due to repeated collisions between 
particles. In addition, we have verified from the profiles 
of the volume fraction and mean velocity that there is no 
clustering instability [15] for the system sizes and volume 
fractions considered here. In order to contrast our results 
to the equilibrium case, simulations are also performed 
with elastic particles in the absence of applied shear at 
the same volume fractions. 
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Fig. 3: (a) obtained from simulations with elastic spheres 
at equilibrium, (symbols: <j) — 0.50 — o, 0.56 — □, 0.59 — o) 
ipit) decays consistent with t"^'''^ for elastic case, (b) V(^) 
obtained from simulations with sheared inelastic spheres. 
for sheared inelastic spheres decays faster than t~^^'^. Lines 
with slope —5/2 (dash dot), —7/2 (dot) and —15/4 (dash) are 
also added to guide the eye. An exponential fit (solid line) is 
also added for comparision. 



Figure 2(a) shows ip{t) defined as, 



{Au,{0)Auj{0)) 



(1) 



where Au is the fluctuating velocity of the particles, and 
i and j correspond to the directions x, y, or z, and the av- 
erage is over all the particles tracked over all times in that 
direction. ip{t) corresponding to the experiments shown 
in Fig. 2(a), is measured by calculating Au from the dis- 
placement of the particle over a time interval over which 
the mean flow travels a distance of O.ld at that location. 
In the experiments, each particle is tracked as long as its 
center is within the layer (Id wide) in the shear direction 
so has to have strain rate approximately constant. We ter- 
minate the trajectory even if the particle leaves the layer 
momentarily. Therefore if it returns back to the same layer 
at a later time then it is considered as a new trajectory. 
ipit) clearly decays rapidly over a time scale less than I/7. 
It is also observed that ipyy decreases more rapidly in the 
boundary regions where particles are layered. ilj{t) would 
collapse on to one curve if T was uniform in the system. 
Because there is an energy current in the experiments, 
and T is not exactly proportional to 7^, multiple curves 
are observed. 

Next, in order to contrast with correlations observed in 



unsheared granular flows, we replot ?/;(i) for the sheared 
case with the corresponding ^(t) for uniform (unsheared) 
granular flow obtained at various flow rates published pre- 
viously in Ref. [5]. Here, we have averaged over the x and z 
directions to obtain smoother curves because ^p{t) appears 
to be isotropic in the vorticity plane. Scaled this way, ip{t) 
collapse on to two sets of curves depending on the pres- 
ence or absence of shear^. In case of the unsheared case, 
ip{t) crosses zero and becomes negative before decaying to 
zero for all the curves corresponding to different flow rates 
(see Fig. 2(b)). Whereas ilj{t) for the sheared case decays 
differently and faster to zero than The same data 

is plotted in log-log scale in Fig. 2(c) to show case the fast 
decays further. It may be noted that the decay is closer 
to the exponential flt plotted than to the unsheared case. 
While fines with slope -15/4, -7/2 [8] and -5/2 [6] are 
drawn to guide the eye, these decays are so fast that the 
range over wfiich ipit) can be measured does not allow us 
to draw any firm conclusions on the scaling of the decay. 

In Fig. 3, ip{t) obtained from simulations of inelastic 
sheared spheres with e„ = 0.9 is shown, along with sim- 
ulations of elastic spheres at the same (j). Decay propor- 
tional to t^"^/^ is clearly observed for low volume fractions 
ip — 0.2, which is not shown in Fig. 3(a) for sake of clar- 
ity. For (j) — 0.4 and higher, it is well known that there 
is a reversal in the sign of tpit) from positive to negative 
at a finite time, due to the reversal in the velocity of the 
particle trapped in the cage between its neighbors. The 
time for this reversal decreases as (j) increases (as can be 
noted from the shift of the negative peak in Fig. 3(a)), 
due to a decrease in the dimensions of the cage formed by 
the neighboring particles. However, the absolute value of 
^(t) still shows a power law decay proportional to t~^^-^ in 
the long time limit. By contrast, the sheared simulation 
data does not show this reversal and shows a fast time 
decay for all (j). Interestingly, it may be noted that ip{t) in 
the three directions are almost equal to each other except 
when layering is important. 

To examine any link between orientational order and 
ip{t), we examine the two dimensional orientational order 
parameter =< exp {i60p) >, where < . > is the average 
over all nearest neighbor bonds and Op is the angle between 
particle bonds. For disordered packing, qq = 0. Evaluat- 
ing qq for the experimental data in the vorticity plane, we 
find that qe decreases from about 0.075 near the side walls 
to less than 0.03, a few grain diameters away from the side 
wall, indicating disordered packing, gg in the shear plane 
is observed to be ~ 0.1 but still significantly smaller than 
qe — 1 for hexagonal order, obtained from the simula- 
tions, both in the vorticity plane and the shear plane, is 
also observed to be less than 0.03, indicating that there is 
no hexagonal ordering. 

However, the state of order appears to be linked to 
particle diffusion in a shear flow of inelastic particles, as 



^It may be more natural to scale with the square root of the 
granular temperature, but at such high densities, it is impossible to 
obtain instantaneous velocities of particles in the experiment. 
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Fig. 4: The scaled mean square displacements (MSD) in (a) the 
flow direction, (b) the vorticity direction, and (c) the shear di- 
rection in the experiments (symbols) and the simulations (solid 
line) for various strain rates as a function of the scaled time f-y. 
The contribution of the affine deformation of the flow has been 
subtracted from the displacements (see text). The dashed line 
corresponds to a slope of 1. 



quantified by the mean square displacement (MSD) of 
the particles (see Fig. 4.) To calculate MSD in a shear 
flow, one has to take into account the affine deforma- 
tion of the flow (sometimes also refereed to as Taylor 
dispersion [16].) Therefore, ((Aa:)^) = {{x[t) -x''{t)f), 
{{Ayf) = m)-yiO)r), and ((Az)^) = ((z(t) - z(0))2), 
where x°'{t), corresponds to the affine deformation of the 
shear flow. After accounting for this contribution, we find 
that slope approaches 1 over long times in the vorticity and 
the flow directions, showing that the diffusion constants 
scale with strain rate. However in the gradient direction, 
MSD in the experiments appear to grow sub-diffusively 
and particles even appeared trapped in a layer over the 
observation window. This effect is most pronounced in the 
region next to the side wall, and the slope only appears to 
approach 1, as one examines regions away from the side 
walls. Because of the Lees-Edwards boundary conditions 
used in the simulations, there is no layering, and a clear 
crossover to diffusive behavior is observed even in the gra- 
dient direction for 'jt/d > 1. 



Therefore, it appears that the rapid decay of ')p{t) is not 
related to the sub-diffusive nature of the MSD (caging ef- 
fect), but appears to be present independent of the state 
of order and the diffusive/sub-diffusive nature of particle 
motion. These results clearly show that the velocity au- 
tocorrelation function decays much faster than the t~^/^ 
decay for a fluid of elastic particles both in the presence 
and absence of ordering, indicating that it seems to be an 
inherent property of a dissipative sheared system. 

Finally, we argue that even though the experiments are 
done with particles immersed in a fluid, the very small ra- 
tio of the viscous and the gravitational forces means than 
the momentum is essentially transfered between particles 
during contact and very little is transferred to the inter- 
stitial fluid just as in dry granular matter colliding in air. 
Therefore, we anticipate that the analysis performed with 
ring-kinetic approach [7] carry over because momentum is 
conserved and energy is dissipated except perhaps in the 
short time limit depending on the details of the interac- 
tion. 

In summary, we have shown that the decay of the veloc- 
ity autocorrelation function in a sheared granular flow is 
qualitatively different from that in an elastic fluid at equi- 
librium. Our results cannot verify the scaling laws [6,7], 
because one would need to track the decay in the autocor- 
relation function over more than two to four orders of mag- 
nitude for one magnitude of increase in the time. However, 
the fast decay in the velocity autocorrelation function in- 
dicates that the dynamics in a sheared granular flow is dif- 
ferent from that in an elastic fluid at equilibrium. Specifi- 
cally, the non-analytic contribution to the shear stress pro- 
portional to 7^^^, which is present in an elastic fluid due 
to long time tails in the autocorrelation functions, may 
not be present in a sheared inelastic fluid. This raises the 
prospect that kinetic theory calculations may be applica- 
ble for a larger range of volume fractions than equivalent 
theories for fluids of elastic particles. Therefore, consti- 
tutive relations based on kinetic theory for granular flows 
may not be restricted to the dilute limit, as previously 
thought, but may be applicable even at higher volume 
fractions. 
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